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Abstract. The problem of separation of variables in some coordinate systems obtained with the use of 
L-transformations is studied. Potentials are shown that allow separation of regular variables in a perturbed 
two-body problem. The potential contains two arbitrary smooth functions. An example of a potential is 
considered allowing explicit solution of the problem in terms of elliptic functions. The cases of bounded and 
unbounded motion are shown. The results of numerical experiments are given. 
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Notation. Everywhere below vectors are regarded as column vectors, and are noted by bold letters. The 
sign T placed over the vector or matrix symbol denotes transposition. A quantity evaluated at the initial 
f^) ' moment of physical or fictitious time is denoted by zero superscript: /(O) = / . 
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H ' 1. Introduction 

The integrable cases of motion equations have great practical value. Their significance is 
determined by the fact that with the help of their solutions one can analyze the motion. In a 
number of cases integrable problems are used to construct intermediate orbits [U |2]. One of 
non-trivial examples of integrated systems is the particle motion in a Newtonian field with 
additional constant acceleration vector. This had been investigated earlier by a number of 
authors [SJIUIS] and applied to analysis of space flights with constant jet acceleration. In 1970 
this problem had been studied using regular coordinates obtained from the KS-matrix [B]. In 
contrast to [6], in [7] integration of the same problem was performed in regular coordinates 
obtained with the use of L-transformations. 
>■ ■ In the present work we consider a problem of constructing potentials allowing integration 

of the equations of motion. The idea of our approach consists in the following. First, a new 
dynamic system is constructed, having more degrees of freedom than the original one. To 
do this, an L-transformation is applied. The theory of L-matrices and their applications is 
given in [HI E] - Using new coordinates, a general potential is selected, allowing separation 
of variables in the Hamilton - Jacobi equation. After this, an inverse transform to original 
coordinates is performed, using explicit formulas. As a basis for selecting general potential 
with the required integrability property, a well known Stackel theorem is used [TU]. This 
^ theorem gives necessary and sufficient conditions for separation of variables for orthogonal 

Hamilton systems, i.e. systems whose Hamiltonian contains only squares of generalized 
momentums. 

Note that separation of variables depends on a choice of a coordinate system. We consider 
here three kinds of coordinate systems: regular, bipolar and spherical. The last two systems 
are introduced in regular coordinates. Canonical equations in regular coordinates are con- 
structed using arbitrary L-transformations from the initial canonical motion equations of the 
perturbed two-body problem. The new equations have also orthogonal form and are invari- 
ant with respect to L-similarity transforms. In the nonperturbed case these equations do not 
have singularity at the attracting center. Due to invariance with respect to some perturbing 
potentials allowing integrability, one can introduce two additional angular parameters. 

As a result of this approach the general solution of original system is represented in 
parametric form, where fictitious time plays the role of parameter, while the physical time 
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depends on this fictitious time and initial data. This sort of integrability is sometimes called 
'Sundman integrability' [TTj . 

As an example of integrable case of the perturbed two-body problem the special kind 
of potential is given. In this example the explicit solution of a problem in terms of elliptic 
functions is expressed, and the criterion of bounded motion is formulated. 

2. The separation of variables 

Consider the Hamiltonian function of the perturbed two-body problem 

H = F(x,y) = ^|y| 2 - - + V, // = 7(771 + m ), r = |x|, (1) 
2 r 

where x = (xi,x 2 ,x 3 ) T is the position vector of the point of mass m with respect to the 
point of mass m ; y = (yi, y 2 , y^) T is the generalized impulses {y,i = ii, i = 1,2,3); 7 is the 
gravitational constant; V = V(x) is the perturbed potential. 

For construction of the equations of motion in regular coordinates we shall need the 
L-transformation z = L(q)q generated by the L-matrix of the fourth order that has the 
following properties: 

L(q)L T (q) = L T (q)L(q) = |q| 2 £ V q G R 4 , (2) 

(L(q)p)i = (L(p)q)i, i = l,...,p, (3) 
0L(q)p)i = -0L(p)q)<, i = p + l,...,4 (4) 
Vq,p G R 4 . 

Here E is the unitary matrix. The conditions (J2J) — (jlj) simultaneously hold only for p — 1 
or p = 3. The quantity p is the rank of L-transformation. The following theorem can be 
proved 

THEOREM 1. An arbitrary L-matrix generating L-transformation of rank three, has 
the form 

( q T K 1 K 4 \ 
q T K 2 K 4 

\ q T ^4 / 

where orthogonal skew- symmetric matrices Ki, K 2 , K 3 , are equal to either 



(5) 



Ki = auU + a 2i V + a 3i W, i = 1,2, 3, , . 

iT 4 = aiAf + a 2 y + a 3 Z, l ° j 

or 

i<Ti = ai»Af + a 2 i3^ + 03^, z = 1, 2, 3, , , 

iT 4 = aiU + a 2 V + a 3 W. [> 

The triplet of vectors = (a^, a 2 j, a3j) T , z = 1, 2, 3, forms an orthonormal basis in R 3 , and 
e = (ai,a 2 ,a 3 ) T is an arbitrary unitary vector. 

Conversely, the arbitrary four skew-symmetric matrices in the form ([6} or <^ define the 
L-matrix by the formula (jSJ). 
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In the formulae 
trices 
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and (E]) there are the so-called basic skew-symmetric orthogonal ma- 
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The matrices are called generators of the L-matrix. If Ki, K%, K 3 , K4 are calculated by 
the formulae flBJ) then L(q) is called the L-matrix of first type, otherwise the L-matrix of 
second type. 

We transfer from variables t, Xi, y% to the new variables r, qj, pj by the formulae 



dt = r dr. 



A(q)q, 



2q 



r A(q)p, q, p G R 4 



(8) 



where the matrix A(q) is found from (jSj) by rejection of the fourth line: 

A(q) = 



q T tfi#4 
q T # 2 #4 
q T K 3 ^ 4 



Consider the equations of motion in new variables qi, Pi 

dlC dp^ _dJC 

dpj ' dr dqj 

with the Hamiltonian 



dr 



-T L = ^~, "P = -t^, j = 0,l,2,3,4 



(9) 



(10) 



K = glPl 2 +Po|q| 2 + |q|V c (q), K(q) = V(x(q)). 

In this system the first equation with j = corresponds to transformation of time: 
|q| 2 c?r. The variable p$ is conjugate to go an d has a constant value. 
If 

x(0) = x°, y(0)=y° 

are initial conditions for the variables of the system with the Hamiltonian (TO , then, as it is 
proved in [T2J Q2] , with the initial values defined by formulae 



?o(0) 
Po(0) 



0. 



A(q°)q°, 
2A r (q°) y ° ) 



(12) 



the solution of (jUJ) becomes, under the transformation ((SJ), a solution of the system with 
the Hamiltonian (jTJ satisfying the initial conditions ( TTTj) . The function q T i^ 4 p preserves a 
constant value along solutions of (JSj), and with the initial conditions from (|12|) . this value is 
zero [13]. Hence, the equality q T i^4p = is the first integral of this system. The variable go 
coincides with physical time t. 
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Note that the systems with Hamiltonian ([I]) and f lTUj) have different orders. The choice of 
initial values by the formulae ( JT2"]) means that there is a special construction of the system 
for each trajectory of the system with Hamiltonian ([T]). 

Let's pick up the form of potential V, admitting division of variables. For this purpose 
we shall take advantage of the theorem proved by Stackel |10j . 
THEOREM 2. The system with Hamiltonian 



n ^ 



admits separation of variables in the Hamilton - Jacobi equation if and only if there is a 
nonspecial matrix <E> of order n wose elements ip S i depend only on q i} such that 



;i,o,...,o) 3 



(13) 



where c = (ci, c 2 , . . . , c n ) T . 

In this case the integrals of motion will be 



t-Pi 



j2 J yii(qi)dqi _p = J2f <Psi(qi)dqi 
y/fi(Qi), i = l,...,n, 



(14) 



Pi 



where = 2(ai<p li (g i ) + . . . + a n (p ni (qi) - Vi(ft)); ft (i = 1, . . . , n) is constant. As g? 

a simple root of the function /i(g.j) is taken. 

Consider again the separation of variables in regular coordinates q^. The Hamiltonian 
looks like f fTUj) . In this case we have 



C\ — C2 — C3 



1 1 4 

c 4 = t, |q| 2 (po + K(q)) = 7 



The solution of system f[T3"j) will be, for example, the matrix 



/ 4 \ 
-110 
0-1 10 
\ 0-11/ 



The potential is defined up to a constant. As po is a constant, we obtain 

1 



K(q) 



4|q| 



;(Vi(gi) + U 2 (g 2 ) + V 3 (q 3 ) + V 4 (q 4 )). 



(15) 



Let's find expression for the potential V c (q) in original coordinates x. Let's notice that 
variables x« and r are quadratic forms of the variables gi, g 2; #3; 94- Using the L-similarity 
transformation it is possible to choose an L-matrix such that a linear combination B1X1 + 
B2X2 + -B3X3 it will be equal to the sum of squares of qi with some coefficients. Note that for 
any L-matrix we have r = |q| 2 . As V c (q) is to be of the form (Tl5i) . the required potential in 
x-coordinates will be the function of the form 



V(x) = -(Ar + Bxxi + B 2 x 2 + B 3 x 3 ). 
r 



(16) 
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Let's specify a choice of L-matrix with the required property. Introduce the notation 



B = yjB* + B* + Bl h = ^, 1 = 1,2,3. 

Suppose that the L-matrix is of the first type. That is, Ki, K 2 , K 3 are calculated by the 
formula ([H]); for simplicity we assume that K 4 = —y. Then 

Ar + B(biXi + b 2 x 2 + 63^3) = Ar - £?q T (biau + b 2 a\ 2 + b 3 a 13 )U+ 



+(6i0 2 i + b 2 a 2 2 + b 3 a 23 )V + (M31 + 62^32 + b 3 a 33 )W 



3>q. 



Choose the parameters ay of L-matrix in such a way that the following equalities hold: 

M11 + b 2 a 12 + 6 3 ai 3 = 1, 

M21 + b 2 a 22 + b 3 a 23 = 0, (17) 
61031 + M32 + 6 3 a 33 = 0. 

Geometrically, the solution to this system means that the vector ii = (an, ai 2 , ai 3 ) T coincides 
with b = (61, b 2 , b 3 ) T , and the vectors i 2 = (a 2 i, a 22 , a 23 ) T , i 3 = (a 3 i, a 32 , a 33 ) T are orthogonal 
to b. Moreover, it follows from the structure of the L-matrix that vectors ii, i 2 , and i 3 form 
a frame. It is evident that the system (|17|) has infinite number of solutions. We write its 
general solution. For the first vector we have 

ii = (b 1 ,b 2 ,b 3 ) T . 

For i 2 and i 3 we assume, in the case b\ + b\ 7^ 0, that 

1 2 = — ^=^= ( 6 2 cos a + bib 3 sin a, —61 cos a + 6 2 6 3 sin a, —{b\ + 6 2 ) sin a ) , 

\Ai + &2 V ' ' 

1 3 = — ^=^= ( — 6 2 sin a + 6i& 3 cos a, 61 sin a + 6 2 6 3 cos a, — {b\ + 6 2 ) cos a ) . 

V b i + h l V " 7 

If b\ + b 2 = 0, then b = (0, 0, 6 3 ) T , 63 = ±1. Therefore, we can take the following vectors as 
the general solution of the system (JTTJ) : 

ii = (0,0,6 3 ) T , i 2 = ^cosa, sin a, , i 3 = 6 3 ^— sin a, cos a, . 

The quantity a G [0, 2ir] plays the role of an arbitrary parameter of the general solution. 

After choosing the parameters a^-, the matrix A(q) is determined uniquely. The solution 
of ( TT7]) gives 



q 



' ((A + B)q\ + (A + B)ql + (A - B)q* + (A - B)fi). 
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|q 

Hamiltonian in g-coordinates corresponding to this potential becomes 



11 

/c = Ei(2^ +4 ^ 2+4A%2) ) ) 

i=l 



where D\ = D 2 = A + B, D% = P 4 = A — B. The canonical system of the equations falls 
into four subsystems 



dr 



TPi 



dpi 
dr 



-2(p + D t ) qi 



1,2,3,4. 



(18) 



These systems are equivalent to four harmonious oscillators. Integrals of motion are obtained 
either from (JHJ), or straightforward from solving (fl8|) . Thus, separation of variables for 
potential (|T6|) is carried out. 

For regular g-coordinates, we introduce a new coordinate system. To preserve the canoni- 
cal form of equations of motion, we use the canonical transformation with generating function 

* = Pxy/Ox cosQ 2 + p 2 \^Qi sinQ 2 + cos Q4 + p^^/Ol sin Q±- 



We obtain 



a = - 
yi dp 1 

qs dp 3 

P _ 5^ 1 



/ Q^sinQ 2 , 



QicosQ 2 , <?2 = &^ ■ 
Q3COSQ4, q 4 = = y/QlsmQ A , 



yp4 

(pi COsQ 2 +P2 SmQ 2 ) 



P-2 

P3 
Pa 



W~2 

TO 
TO 



y/Qx{—pi sin<5 2 + P2 cosQ 2 ). 
o ]rr- (Ps cos ^4 + P4 sin Q 4 ) , 
v/Q^(-p 3 sin Q 4 + p 4 cosQ 4 ). 



(19) 



The coordinates Qi, Q 2 , Q3, Q4, obtained from (fl9|) . will be called bipolar. From the last 
four equations we find p±, p 2 , P3, p 4 : 



Pi = 2P iy /Q^cosQ 2 



j^=smQ 2 , p 2 = 2P l ^/Q~ 1 smQ 2 + 



cos Q 



2: 



P3 = 2P 3V / Q^cos(5 4 ^ sin Q 4 , p 4 = 2P 3V / Q 3 "sin Q 4 

V V3 



cos Q 4 . 



(20) 



In the new variables the Hamiltonian /C becomes 



+ pMi + Q 3 ) + (Qi + Q 3 )V, 



where function V is expressed in terms of Qi. 

Similar to the above, consider separation of variables in bipolar coordinates. In the 
notations of theorem 2 we now have 



c i — Qii 



c 2 



4Qi 



C 3 = Q3 



1 



4Q 2 



As a solution to (fT5|) one can take the matrix 



¥ 
w 1 





\ 








4QI V 



(21) 



For the potential V admitting separation of variables, we find 

In g-coordinates we obtain the form 

1 / „ „ „ „ V^farctan— ) „ „ „ „ V^farctan— x 

v. - w (M + 4)v^ + d) + + (d + d) v,M + d) + - 1 | IT ^f 

Passing to x-coordinates, we use the concrete L-transformation 

= 2q 1 q i + 2q 2 q 3 , 

-- -2q iq3 + 2q 2 q±, (22) 
= Ql + <£ - <& - Ql, 

which follows from ([5]), (jSJ) with K\ = V, K 2 = W, =U, K4 = —y. Taking into account 
that for any L-matrix the equality r = qf + q% + q 3 + q% holds, we obtain 

1 1 
2,2 L 1 \ 2,2 / 1 \ 

% + <? 4 = ^ r ~ X3 >i Qi+Q2 = 2^ r + X3 ^ 
The general solution of the first equation is 




<?3 = \j r ^ cos ^' 94 = y r 2 Xs sin?/', V e [0, 27t]. 
Then 

x\ sin ip — x 2 cos "0 x\ cos ?/> + 22 sin ip 

V2y/r - x 3 " v^ 7 " - ^3 

In a similar way we may introduce a parameter, using the second equation, 



qi = y r cos^i, <?2 = y r sin^i, ^1 e [0, 2tt]. 

As is well known [12] , under L-transformation for a point in R 3 at a distance r from origin, 
there corresponds a point of some circle of radius ^Jr in R 4 . The variables q^ contain an 
arbitrary parameter ip (or ^1), giving parametrization of the given circle. In the original 
coordinates Xi this parameter disappears. Note that 

— = tan ipi, — = tanip. 
Qi 13 

We therefore assume functions V 2 , V4 to be constant. Then we arrive at a potential of the 
form 

V(pc) = - k((r + x 3 )/2) + G 2 ((r - x 3 )/2)} , (23) 
r i 

where Gi, G 2 are arbitrary smooth functions. The Hamiltonian in bipolar coordinates for 
this potential takes the form 

/C = Q 1 ( T +Po + ^^j + i ^ T + Q3( T + P o + ^ r j + i ^ T . 



In view of the solution ( 12 ip for /j from the theorem 2 we have 

G 1 (Q 1 ) 



Qi 4Qf Q 



«3 

77" -Po 



Po 



Qi 
G 2 (Qs 



h{Q2 



2a 



2; 



, U(Qi) = 2a 4 . 



Then integrals of motion are obtained by formulas (JED- 

Let's consider one more case of separation of variables. Introduce in g-coordinates the 
spherical coordinates 



fQx cos Q 2 cos Qi , q 2 = VQi sin Q 2 cos Q A , 
fQl cos Q3 sin Q4, <?4 = \/£?i sin Q3 sin Q 4 . 



(24) 



We supplement the transformation (|24p to obtain a canonical transformation of impulses 



o /7T~ /O /O d SU1Q2 D cos Q 2 sin Q 4 D 
Pl = 2 VQi COS Q 2 COS Q4P1 /7T~ ~ >, ^2 /7T~ - A 

'Q1COSQ4 vQi 



P2 = 2v ^sinQ 2 cosQ 4 Pi + ™ Q * P 2 ~ **Q'™Q* R 

/QicosQ 4 yQi 



4; 
4, 



P3 = 2v^cos Q 3 sin Q4P1 - ^9* P 3 + C ° Sg ^ sQ4 P, 

VQi sin Q 4 VQi 

P4 = 2yQrsinQ 3 sinQ 4 Pi + ™9* P 3 + "^ff^ 4 ^. 

VQismQ 4 VQi 



Then in new variables the Hamiltonian will be 



1 / P 2 



4QiP^ + 



+ 



p2 
^3 



+ 



p2 
-"4 



Qx cos 2 Q 4 ' Qi sin 2 Q 4 ' Qi 

In the notations of Stackel theorem we have 

1 1 



+ P0Q1 +Q1V. 



ci = Qx, c 2 



^3 



4Qi cos 2 Q4 ' ~° 4Qi sin 2 Q 4 ' 
In this case the solution of (fT3l) will be the matrix 



c 4 



4Qi 



4 

1 

1 





\ 
1 



cos 2 Q 4 
1 

sin 2 Q 4 



The potential V, admitting separation of variables, can be written as 

V = ^-(QiV^Qx) + 1 F 2 (g 2 ) + 1 F 3 (Q 3 ) + t^-v 7 4(Q4; 

V AQi cos <y 4 4Qi sm z Q4 4Qi 

In view of relations 



(25) 



(26) 



/o 2 -o 2,2 r + ^3 /o • 2 /n 2,2 r x 3 rs 11 

Qi cos Q 4 = 9i + q 2 = — - — , Qi sin Q 4 = q 3 + <? 4 = — 7, — > Qi = IqI 



i 2 Q4 = % % = —r-i Q2 = arctan — , Q 3 = arctan — 

q( + q 2 l + x 3 /r ' q 1 q 3 



following from ( |22|) . ( 124"|) . and the remarks above, we obtain the required form of potential 
in x-coordinates 



V{x) = - G 1 (r) + ^- + ^- + -G 2 (^)], (27) 
r + x 3 r — x 3 r r J 



1 

r 

where Gi, G2 are arbitrary smooth functions and A, B arbitrary constants. 

Now assume that a Hamiltonian (pQ) with the potential (|27|) is given. Applying L- 
transformation (|22p . we write the new Hamiltonian in g-coordinates as 

^ 1| 12 1 12 ^ /: :2\ ^ £? 1 _ (q\ + gf - g 2 - q\ 

8 qf + qi qi + qi |q| 2 

Fulfilling canonical transformation (I24p . (125]) . we have 



2 ^ u Qi / 4Q X cos 2 Q 4 V 2 

+ ^ 1 27; f^f + 45 ) + 77T + 4G 2(cos 2Q 4 ; 
4^! sir Q 4 V 2 / 4Qi V 2 

Taking into consideration matrix (l26p. we then obtain 

f 3 (Q 3 ) = 2(« 3 - 45), / 4 (Q 4 ) = 2 ( - a 4 - 4G 2 (cos 2Q 4 )) . 

\ COS Lg4 Sill t^/ 4 / 

The integrals of motion follow from (fl4l . 

Note that using arbitrary L-transformations allows to introduce two parameters into the 
potentials obtained. Tthese two parameters are determined by some constant unit vector b. 
For example, instead of (|2"T|) one can write 



K(x) = - 

r 



■ , , 2A 25 1 /b J x 

Gl r + jTTTt + ij" + -G 2 ( 

r+b J x r — b J x r Vr 



In the next section we show how to perform separation of variables in this case. 

3. Integration of the system of equations in a special case 

In this section we perform straightforward integration of a system with potential of the 
form (123]) having additional parameters. Namely, consider the potential 

V = V(x) = - l - (g x ({t + b T x)/2) + G 2 ((r - b r x)/2)) , (28) 

where Gi, G 2 are some smooth functions, and b = (pi, b 2 , b 3 ) T an arbitrary unit vector. Note 
that the vector b provides two parameters in explicit form. Having in mind only theoretical 
investigation (integrability problem), one can take b to be the ort along the xi-axis. On the 
other hand, from the more practical point of view, introducing vector b gives us additional 
degree of freedom necessary for applied problems of celestial mechanics. In such problems, 
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the axes are usually connected with some special directions (equinox or zenith). Therefore 
the presence of the vector b in potential (|28p allows one to turn the coordinate system at 
one's will. 

As G±, G 2 one can take, for example, functions of the form 

V + b T x) fc , -V-b T x) fc , A; = 1,2,... 

We consider a finite linear combination 

1 N 

V=—J2 {Mr + b T x) fc + B k (r - b T x) fe ) . (29) 

k=l 

Here A k , B k are constants. Such a potential was considered in [16]. This case leads in general 
to hyperelliptic integrals. 

For an interested reader here is a problem: find a real perturbing potential which can be 
approximated by functions of the form f[29~l) . Note that the combination 

-E( r + b T x) 2 + — (r - b T x) 2 = -5b T x 

gives potential corresponding to a constant force. Applications of such potential were con- 
sidered in [3J IH IS] • 

The canonical equations of motion have the form 

it = Vh 

|?= -^-^(G 1 ((r + b r x)/2) + G 2 ((r-b^x)/2))+ (30) 

+ 27 ( G 'i(( r + bTx )/ 2 )(f + h i) + G 2( r ~ b T x)/2)(f - 6,)) , 

where i = 1,2,3 and the sign prime indicates the derivative. 

This system is the same as the equation of the perturbed two-body problem 

* + F X = i ( G 'i ((r + bTx)/2) " G ' 2((r " bTx )/ 2 )) b + 

+ 2^2 (<?i((r + b T x)/2) + G' 2 ((r - b r x)/2))x- 

(G x ((r + b r x)/2) + G a ((r - b r x)/2))x. 

From this one can see that the perturbation is defined by two forces. The first force is 
collinear to the fixed vector b, and its module varies in dependence on vector x. The second 
force is the central one. 

We are going to show that the system (I3"0"|) is integrable in regular variables found by 
L-transformations. Transformation (JS]) contains an arbitrary L-matrix. A special choice of 
this matrix allows one to separate the variables in the case of an arbitrary constant unitary 
vector b. 

Consider the term in fflU]) containing V c (q). In the new variables this becomes 
|q| V c (q) = "^((Iql 2 + q T (6 1 ^ 1 + b 2 K 2 + b 3 K 3 )K,q) /2)- 
-G 2 ((|q| 2 - q T (&iKi + b 2 K 2 + b 3 K 3 )K 4 q)/2). 
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We assume that the L- matrix has the first type and K4 = — y. Then 



|q| V c (q) = -G x {{\q\ 2 - C)/2) - G 2 ((|q| 2 + C)/2), 



(31) 



where 



(Mil + M12 + b 3 a 13 )U+ 



C = q T 

+(b 1 a 2 i + b 2 a 2 2 + b 3 a 23 )V + (&ia 3 i + b 2 a 32 + b 3 a 33 )W 
Let's select parameters L-matrixes from a system (j!7p . Then 

C = q T Uyq = q T 



yq. 



-92 

93 



2 2,2,2 

5l -52 + ?3 + 5 4 - 



V 54/ 



Substituting the found value C in (T5TT) . we obtain 

|q| 2 K(q) = + q\) - G 2 (q 2 3 + q 2 ). 

It follows that the Hamiltonian ffTUj) is represented in the form of the sum 

K = fC\ + /C2, 

where 

£1 = g (P? + Pi) + Po(ql + q\) ~ G x {q\ + q 2 2 ), 

£2 = l(pi + pI) + Po(q 2 3 + 5 4 2 ) - ^2(53 + 5 4 2 )- 

As the value of po is constant, the system (jHJ) splits into two independent subsystems 

dqi dfCi dpi dfCi 



dr 



dpi 



dqi _ dJC 2 
dr dpi 



dr 

dpi 
dr 



dqi ' 

dJC 2 
dqi ' 



z = l,2, 



3,4. 



(32) 



(33) 



We integrate all over again a system (1321) . In the bipolar coordinates Hamiltonian K.\, 
and accordingly the system, have the form 

K x = \{±QiPl + ^) +P0Q1 - Gi(Qi), 



dr 



/n p <^t?2 _ P 2 

1 P 2 

1 D2 ^ 



dP 2 



+ gg2 _ Po + G[(Qi), -j^f 



0. 



(34) 



Since the Hamiltonian /Q does not explicitly depend on r and Q 2 , the system (134)) has 
two integrals, 

\QiPl + ^ + P0Q1 - GxCQi) = y • (35) 
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^2 = Ci. 

Here, Ex and c\ are the constants of integration. Taking these integrals into account, the 
equation for Pi may be written in the following form 

dr ~ AQ\ 8Q1 liQl) Qx ' 
Eliminating dr from equations for Pi, Qi and integrating the resulting equation, we find 

Si 



where 

$i(Qi) = -c\ + ExQx + c 2 Ql + 8QiGi(Qi) 
and c 2 is integration constant defined by 



c 2 = 4(P{ 



0,2 C l ^1 a Gl(Ql) 



(QI) 2 Qi Qi ' 

Due to nonnegativity of Qi, from the first equation of the system (|34|) it follows that 

5i = signQi- 

Substituting the derived Pi to the first equation of we find 

Qi 

T + c a = 28 1 f-^==. (36) 

Using the continuity principle, the sign before the integral (|36|) cannot change when $i(Qi) 
is non-zero. Therefore, the function t(Qx) in this case behaves monotonically. Inverting 
the integral (1561) . we obtain Qi as a function of r; we substitute this function in the second 
equation of the system (13"4l) . Then we get 

T 

Ci f dr q 



o 

Here c 3 and c 4 are the integration constants. Thus, the values Qx, Q 2 , Pi are represented as 
functions of r. If <&x{Qi) is a polynomial, the integral (13"6l) is, in general, hyperelliptic. 
The integration of the system (1531) is done similarly. We find as a result 



Qi — ~r I r, , x + Cs, c 8 — <5 4 , 


^3 = 7^- VMQz), $2 = sign Q' 3 , Px = c 5 , 

where 

MQs) = -4 + E2Q3 + c 6 Ql + 8Q 3 G 2 (Q 3 ). 
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Here c§ and P 2 are the integration constants defined by the equalities 



c 6 = 4(P 3 °) 2 + 4 E2 ° GM) 



(Q3) 2 Q3 Qs 



E > 1 Q 3 P 2 + ^r+ P0Q3 - G 2 (Q 3 ). (37) 



T " 2°*^ + 8Q 3 
The function ^(r) is found by a reversion of the integral 



r + c 7 = 25 2 /^£%=. (38) 
n 

Thus, the values Q 3 , Q 4 , and P3 are also determined as functions of the variable r. The lower 
limits £ and 77 in integrals ([3T)j) and f[3"8l are chosen according to the location of Qi and Q 3 
with respect to the roots of functions <J>i(<5i) and ^ 2 (Q 3 ), respectively 
The formulae of inverse transformation 



Qi = qf + ql tanQ 2 = | 2 -, 
Q3 = qi + qh tanQ 4 = ||, 

allow to define initial values of the variables Qi and Pj (i = 1, 2, 3, 4). 

The values of integration constants c\, c 2 , c 3 , e 4 , and Pi are determined by the initial 
values of Qi, Q2, -Pi, P2- These five constant values are connected with each other by the 
integral f l35]) . In the same way, the constant values c 5 , c 6 , c 7 , c 8 , and P 2 are connected by the 
integral fl3T|) and are defined by initial values of Q 3 , Q 4 , P 3 , and P 4 . From p = — P(x°,y°) 
we find also relation Pi + E 2 = S/j,. One has to add the above relations for c 2 and Cq to these 
connections. Besides, as the bilinear relation q T P 4 p = is the integral of ([9]), in our case we 
have 

q T (-F)p = -q 2 pi + qip 2 + q*p 3 ~ M>4 = 0. 

Therefore the equality P 2 = P 4 , or equivalently c\ = c 5 , also holds. 

Applying further the first four formulas (|T9j) and (1201) . we find g», pi (i = 1,2,3,4) as 
functions of r. Finally, integrating the two remaining equations of ([9]), we obtain po = 
— P(x°,y°) and physical time expressed through r, 



t = qo = I |q| dr + c 9 = ti + t 2 , (39) 

where 







h = J Qi(r)dr, t 2 = J Q 3 (r)dT, c 9 = 0. 


Thus, the system ([9]) is completely integrated and we can, at least locally, find a required 
trajectory. Here it is necessary to note, that if perturbing potentials G\, G 2 in (1301) are 
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analytic, then, as it is known from a course of the differential equations, the solution of a 
problem will also be analytic. Let us suppose that the local inversion of integrals (ESJ, (J3"8]) 
appeared to be a globally determined function. In this case we can conclude, by uniqueness 
of analytic continuation, that this inversion gives not only local, but also global solution of 
the problem (130]) . This is the case when functions G\, G 2 are polynomials of degree two or 
three. In this case ( 136]) and (138]) are the elliptic integrals, for whose inversion we have at our 
disposal the well developed technique of elliptic functions; thus, we have found the solution 
of ( 13 Op in explicit form. 

4. Inversion of the integral in elliptic case 

In this section we consider one case of functions G\ and G2, which reduces to elliptic integrals. 
Other cases have been studied in [TJ HU [161 EH] • Take as G\ and G2 the functions 

Gi= A Tl +A 1 (r + b r x) + A 2 (r + b r x) 2 , (40) 
r + b-' x 

G 2 = B ~ 1 T + E 1 (r-b r x) + E 2 (r-b r x) 2 , (41) 
r — b 1 x 

where A±, A%, B-i, Bi, and B2 are the parameters of the potential. Then for the 

functions $i(Qi) and ^2^3) in ([36]) and (138]) we have the expressions 

$i(Qi) = ci + + c 2 Q 2 + 32A 2 Q?, 

$ 2 (Q 3 ) = c 5 + £ 2 Q 3 + c 6 Q 2 3 + 32B 2 Ql 

Here ci = -c\ + 4A_i, c 2 = \§A X - 8p , c 5 = -c\ + 4S_ X> c 6 = \§B X - 8p . 

Firstly, let us note that the variables Qi and Q3 are non-negative by definition, and that 
from integrals (1361) and (I38p it follows that the ranges of these variables are determined by 
the inequalities 

$i(Qi)>0, $ 2 (Q 3 )>0. (42) 

Let us reverse the integral (136|) . The number of roots of the polynomial $1 and their positions 
depend on the value of A 2 . With A 2 = the degree of $i(Qi) equals to two. The integral (1561) 
is found in elementary functions, so this case is not being considered here. We distinguish 
two cases: A 2 < 0, A 2 > 0. Let's note the roots of as £1, £ 2 , £3. The cases under 

consideration will be sequentially numbered by parameter %a- 

I. Assume that A 2 < 0. In this case 00) > 0, $i(+oo) < 0. The value $i(0) = ci 
may be both positive and negative. For actual motion there should be at least one positive 
root. The qualitatively different cases of the graph of $i(Qi) are shown in figures [T] and [2J 
In the case of three real roots (fig. EJ), the axis of ordinates goes between £1, £ 2 if ci < 0, and 
left with respect to £1 or between £ 2 , £3 if ci > 0. 

The case %a = 1- Suppose that $1 has one real root £1, and that Qi G (0,£i) (fig. H]). 
Let's write the integral (I36p in the form 

Qi 

5\ f dz 

r + c 3 



I^IM} y/(^-z)(z 2 +bz + c)' 



Si 
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Figure 1: The graph of $i(<5i)- The case A 2 < 0. 



where the square trinomial z 2 + bz + c has no real roots and is positive for all z, and 

Apply in the integral the substitution 

1 — cos (p 



and put the notations 



1 + COS if 



^ = 2arctanJ^-^, A; 2 = - (l + Y / x = 2 V /Z 2^4 

V a 2 V a / 



1/, , £i + &/2^ 




£2 \£a Qi 



Figure 2: The graph of $i(Qi)- The case A 2 < 0. 



Then we derive 



5l 

r + c 3 = - — 



^i J y/l- kj sin 2 



Putting here r = 0, we find an integration constant C3: 

v? 

sign Pf 1 f dip 



c 3 



^i 7 yjl - k\ sm 2 y? 



, y?? = 2 arctan 



6 -OS 
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Check that k\ < 1. As z 2 + bz + c has no real roots, we have b 2 — 4c < 0. Therefore, 



(^ 1 + ^) 2 <tf + b£i + c = a 2 



6 + 6/2 



< 1. 



Hence, |fei| < 1. Reversing the integral (jHJ) derived above, we come to the function 

2a 



Qi = £i + a 



1 + cn(Zx(r + c 3 ); fci) ' 



It is easily to see that for Q\ 6 (0, £i) the denominator cn(ix) + 1^0. Calculating the 
derivative of Qi, we get 8\ = — signsn (h(r + c 3 ); fci). For the variable Q 2 we find 



Cir aci 
V2 — T77 — ; — c + 



4(6 + a) 2^ -a 2 ) 



«i(t+c 3 ) 



m = l + 



l + Tix cn(u; 
2a 



1 + rii cn(u; ki) 



6 - a 



Note that 



n 



I 7.2 C 



*1 



> 0. 



- 1 1 4a£i 

Therefore for calculating the integral of (1 + ni cn(w; fci)) -1 we apply the formula (341.03) 



l + ncn(u;k) 1 — n 2 L \ ' n 2 — 1 



ft 



(45) 



with 



\/fc 2 + /c /2 n 2 sn(w; fc) + -\/n 2 — 1 dn(w; fc) 



\/k 2 + fc' 2 n 2 sn(w; fc) — \/n 2 — 1 dn(tt; A;) 



A; = 1 — fe . (46) 



Here we note the elliptic integral of the third kind as 

u 

dv 



n(w, n; k) 



1 — nsn 2 (v; k) 



For t\ we have 



*i = (6 + a)r - — 



«l(r+c 3 ) 



du 



llC3 



1 + cn(u; ki) J 1 + cn(w; fci 



d-u 



The integral of (1 + cn(w; k\)) 1 is calculated by the formula (341.53) [15 



dv 

1 ± cn(v; fc) 



u - £(u) ± 



dn(w; fc) sn(w; A;) 
1 ± cn(w; fc) 



(47) 
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where E(u) = E(ip; k) is incomplete elliptic integral of the second kind (<p = amu). 

The case %a = 2. Suppose that has three real roots < 6 < £ 2 < 6> and 

Q° G (0,6)- Let's write (JMD as 



t + c 3 



61 



dz 



Making the substitution (p = arcsin a/ (£1 — <z)/(6 — z ) an d reversing the resulting integral, 
we find 

n p 

Vl — ?2 



cn 2 (/i(r + c 3 ); /ci)' 



where 



c 3 



sign P® 



' 6~6 

dip 



h = V-2M&-Zi) 



(p[ = arcsin 



' 6 - gg 
6-0?" 



1 — A; 2 sin 2 </? ' 

u 

Now we calculate <5i. We differentiate Qi and use the formula of double argument for elliptic 
sine. We have 

Qi = 2Zi(£ 3 - 6)cn" 3 (n; fci)(-l) sn(u; fci)dn(«; fci) = 
= — ii(^3 - 6)cn~ 4 (w; fei)(l - /c 2 sn 4 (w; fei))sn(2u; fei), 
where the notation u = l\{r + c 3 ) is introduced for brevity. Therefore, 

5i = signQi = - sign sn (2/i (r + c 3 ); fc x ). 

Now we find Q2 

6 



= CiT Ci(6 ~6) 

^ 2 46 4/^6 



n(/i(r + c 3 ),rai;/ci) - n(Zic 3 , riijfci) + Q 2 > "i 



6 



For the value of physical time, corresponding to the variable Q\, we have 



h = &r + 



6-6 

h 



h{r+c z ) 



du 



I1C3 



cn 2 (w; hi) J cn 2 (w; k\ 



du 



where the integral from cn 2 (u; k\) is calculated by the formula (313.02) [15 



dv 



cn 2 (t>; k) 1 — k 2 



1 \l-k*)u-E(u) + Mu]k)sn{u]k) 



cn(w; k) 



The case %a = 3. The polynomial $i(Qi) has three real roots 6 < 6 < 6> an d Qi 
(max{0,6},6)- Write (ESJ) as 



T + C 3 



6l 



dz 
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The reduction of this integral to the standard form fj44|) is carried out by the substitution 
(f = arcsin a/ (£3 — z)j (£3 — £2). The result of reversion can be presented in the form 



where the following notations are used 



c 3 



h = 

sign P® 



?3 - 4"1 



dip 



<Pi = arcsm^/- — 

S3 - 4"2 



1 — fef sin 



For 5i we find 61 = — signsn (2Zi(r + C3); fci). Substitute Qi in the formulae for Q2 and £1. 
We find 



Q2 



n(Zi(r + c 3 ), m; fci) - n(/ic 3 , m; fci 



+ g 



2- 



f 1 = 4V + 



/1 



«l(r+c 3 ) 



I1C3 



sn (w; k\)du — J sn (w; k\)du 

The integral from squared elliptic sine is calculated by the formula [15] 



sn 2 (f ; A;)<it> = -r^(u — E(u)). 
k z 



II. Assume further that A 2 > 0. Now we have $i(-oo) < 0, $i(+oo) > 0, and $i(0) = ci. 
The qualitatively different cases of the graph $i(Qi) are shown in figures [3] and |U 





if 1 







Figure 3: The graph Case A 2 > 0. 

The case %a = 4. The polynomial $i(Qi) has one real root £1 and, accordingly Qi(0) > 
max{0,£i}. The graph of $i((Ji) in this case is shown in fig. |3j Write the integral (136]) as 

Qi 

5i f dz 

T + C3 ~2VM = 2 J v/(z-6)(z 2 + bz~+c) 1 
€1 
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Figure 4: The graph Case A 2 > 0. 

where the square trinomial z 2 + bz + c > for all z. The coefficients b and c are found by 
the formulae fj43l) . Applying the substitution 

l — cos (p rz ~ 

z = ^i + a— , a=W£f + &£i + c 

1 + COS if v 

and reversing the resulting integral, we come to the function 

2a 



Qi = £i - a + 



where 



1 + cn(7i(r + cs); fci)' 

*} = I(i-&±!/?), (l = 2v ^. 



c 3 



sign P® 



dip 



1 — k 2 sin 



(^5 = 2 arctan 



QS-6 



As above, one can show that fc^ < 1. The resulting function Qi(r) is unbounded, as it has 
an infinite number of poles on real straight line, which are found by the formula 



Am + 2 

r = — K(ki) — c 3 , m G Z. 



Further we find that 5i = signsn (h(r + C3); k\). For variable Q2 we have 

r-h(r+c 3 ) l lC3 
C-yT fi.r.-i f fJ.il. 

Q2 



4(& -a) 21^1-a?) 



1 + ni cn(u; k\) 
2a 



du 



1 + n\ cn(w; k\) 



Note that 



rr. 



nl-1 



1 + 



ni-1 



61 + a 

(6 ~ a) 2 
4a6 



< < jfcj. 
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Therefore for calculating the integral from the function (1 + m cn(w; fci)) 1 the formula (T45 
is to be applied with 



9i 



1-n 2 



k 2 + k' 2 n 2 



arctan 



k 2 + k' 2 n 2 sn(w; k) 
1 — n 2 dn(w; k) 



k' 2 = 1 - k' z . 



If = then rii = —1, and for calculating Q 2 the formula ( |47|) should be used. For ti we 
have 

■/l(T+C 3 ) hc 3 



2a 
77 



1 + cn(w; k\) J 1 + cn(w; fci) 



Suppose that <£>i has three real roots £i < £2 < £3- The graph of the function $i(Qi) in 
this case is given in fig. HI This case also splits into two subcases: ^1 < Q\ < £2 and £3 < Q®. 
The case i& = 5. Suppose that G (max{0, £1}, 6)- We write (15BT) as 



T + C 3 



6l 



dz 



6 



We apply to this integral the substitution tp = arcsin y(z — 6)/ (£2 — 6) and use the nota- 
tions 

V S3 _ Si 

Then our integral has the standard form 



5l 

T + C 3 = T" 



dy> 



^1 7 \f\-k\ sin 2 



(48) 



Reversing (j4"8"|) and using the inverse substitution, we find the required function 



Qi = 6 + (6-6)sn 2 (/i(r + c 3 );A; 



where 



c 3 



sign P° 



v? 



dip 



'I - k\ sin 2 <y9 



ip{ = arcsin 



6-6 



As above, one can show that Si = signsn (2Zi(r + c 3 ); For Q2 we find 



Q2 



4/i6 



n(/i(r + c 3 ),rai;/ci) -n(/ic 3 ,ra 1 ;/ci) + Q° 2 , m = 1 



For the first summand of physical time i in f[5^|) we have 



h = 6t + 



6-6 



il(r+C3) 



Z1C3 



sn 2 («; ki)du — / sn 2 (w; k\)du 
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The case ia = 6. Suppose Q\ G (max{0,£ 3 }, 00). The integral fl36|) has the form 

Qi 

Si 



r + c 3 



dz 

2V2A2J ^(z-^)(z-^ 2 )(z-^y 



Transforming this integral to the standard form (|48|) is made using the substitution (p 
arcsin y(z~— £ 3 ) j(z — £2). The resulting reversion of the integral in this case is following 



Qi = £2 + 



£3-6 



cn 2 (Zi(r + c 3 ); £4)' 



where 



c 3 



sign P[ 



dip 



v 



n ■ / Qi — £3 
<y9 x = arcsin , 



sin <y? 



OS -6' 



Now the function Qi(r) has an infinite number of poles of the second order, hence it is 
unbounded. The poles are found by the formula 

2m + 1 

t = — K(ki) — c 3 , m G Z. 

'1 

Further we find the values Si, Q2, t\ 

S 1 = sign sn (21 1 (r + c 3 ); fci), 

6 



Q 2 



err ci(£ 2 
46 4/^26 



n(/i(r + c 3 ),ra 1 ;A;^ 

ii(r+C3) 



n(Z 1 c 3 ,ra 1 ; ki 



Tlx 



6' 



h = 6t + 



/1 



ilC3 



cn 2 (u; ki) J cv?(u;ki) 




An inversion of the integral (|38p is fulfilled in a similar way. This integral and the 
function $2 differ only by notations from the integral ( 136|) and the function $1. Therefore, 
after some evident renaming, we find the expressions for Q 3 , Q4, and t^. We number these 
cases sequentially by the parameter %b- Then we have: 

772,773 G C, < Q\ < rji (Q 3 is bounded). 
< Q® < 7]i < I]? < r] 3 (Q 3 is bounded). 
r]i < T] 2 < Ql < rj 3 (Q 3 is bounded). 
772, 773 G C, 771 < Q° (Q 3 is unbounded). 
T)i < Ql < T] 2 < r)3 (Q 3 is bounded). 
rji < r\2 < r] 3 < Q® (Q 3 is unbounded). 
The study above yield the following theorem. 

THEOREM 3. The motion of the particle is bounded if and only if at the initial moment 
both variables Qi and Q 3 are restricted on the right by the roots of the polynomials $1 and $2? 
correspondingly. 



>B 


= 1 




B 2 


< 


IB 


= 2 




B 2 


< 


'B 


= 3 




B 2 


< 


IB 


= 4 




B 2 


> 


IB 


= 5 




B 2 


> 


'B 


= 6 




B 2 


> 
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Now we give a definition of retaining potential, introduced in |14j . 

DEFINITION. A potential is named as retaining, if for arbitrary initial conditions the 
motion of a particle in a perturbed field corresponding to this potential is bounded. 

Thus, potential (j28|) . where G±, G 2 are defined by the formulae (HUi) . (HTi) for A 2 < 
and B 2 < 0, is retaining. Generally, the formulae (136)) . (138]) are not elliptic integrals, and we 
cannot present a solution in explicit form. Nevertheless, the above-stated qualitative result 
remains true [15] . 

5. Numerical examples and analysis of motions 

In the examples below we consider the motion of a particle in perturbed gravitational 
field of a planet with spherical density distribution, whose gravitational parameter is taken 
to be fi = 398601.3 km 3 /s 2 . The perturbing force is defined by the potential (1281) . with G\ 
and G 2 calculated by the formulae (j4"01 and (fldlh For convenience (to have no fractions), a 
dimensionless direction vector b for the constant force is used. While doing calculations, this 
direction vector is assumed to be normalized. The dimensions of parameters v4_i and B-\ 
are [km 4 /s 2 ], A 1 and Bi are [km 2 /s 2 ], A 2 and B 2 are [km/s 2 ]. Calculations and construction 
of orbits were performed using the Maple system with 32 digits. In each example, for con- 
venience of its analysis, the values of circular and parabolic velocities v C i r , v par of Keplerian 
motion are given. The perturbations being considered are great, they are non- typical for the 
Earth's satellites. For this reason, we do not give Keplerian elements of osculating orbits for 
the corresponding initial values. The initial position of a particle is marked by a point on 
the corresponding figure. 

Example 1. Initial values of coordinates and velocities of a particle: 

X\ = 8200 km, x 2 = km, x 3 = 6000 km, 

x 2 = 8.6 km/s, Xi = x 3 = km/s (v cir ~ 6.26 km/s, v par ~ 8.86 km/s). 

In an unperturbed case these values define an elliptic motion. 
Parameters of potential follows: 

A_i = 0.004 km 4 /s 2 , Ai = 0.06 km 2 /s 2 , A 2 = 0.2 • 10~ 7 km/s 2 , 

B_ v = 0.0001 km 4 /s 2 , B x = 0.008 km 2 /s 2 , B 2 = -0.3 • 10~ 4 km/s 2 . 

Coordinates of direction vector are b = (—1, 2, 1) T . In the case under consideration the roots 
of polynomials $1 and $2 are 

6- 1478, £2 - 115346, q° ~ 463 i Q? e 

772 w 1707, r) 3 w 31031, Q° 3 w 5529 Q° 3 e (r/ 2 , rj 3 ). 

Therefore, the motion is bounded. This is the case %a — 5, %b — 3. 

The coordinates and velocities have been calculated during a time range, corresponding 
to two revolutions of the particle around the attracting centre without perturbations, that is 
r G [0, 2T], where T is calculated by the formula 
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100000 




Figure 5: The case i& = 5, ig = 3. 

Here hf. is the Keplerian energy. Let's remind that L-transformation doubles the angles at 
the origin of coordinates. 

Note that in this example the potential is not retaining. Nevertheless, the motion appears 
to be bounded. The trajectory of the particle is shown in fig. 

Example 2. Initial values of coordinates and velocities of a particle: 

xi = 8200 km, x 2 = km, x 3 = 6000 km, 

x 2 = 9.9 km/s, Xi = x 3 = km/s {y cir ~ 6.26 km/s, v par ~ 8.86 km/s). 

In unperturbed case the motion belongs to hyperbolic type. 
Parameters of potential follows: 

A_x = 0.004 km 4 /s 2 , Ax = 0.006 km 2 /s 2 , A 2 = -0.2 • 10" 7 km/s 2 , 

S_! = 0.0001 km 4 /s 2 , B 1 = 0.008 km 2 /s 2 , 5 2 = -0.3 • 10~ 7 km/s 2 . 

As A 2 and _B 2 are negative we have a retaining potential. Coordinates of direction vector are 
b = (1,2, — 1) T . The roots of polynomials: 

£2 « 2126, £3 « 122192633, Q? « 5529 =^ Q? G (6, 6), 

^ 2 « 1699, r/ 3 « 81506371, Q° w 4631 Q° 3 e (r/ 2 , r/ 3 ), 

Therefore, the motion is bounded. This is the case — 3, ig = 3. The integration is carried 
out during the time range corresponding approximately to t — 1759.74 days. The particle 
trajectory is shown in fig. [6l 

Example 3. Initial values of coordinates and velocities of a particle are as follows: 

X\ = 6000 km, x 2 = km, x 3 = —8000 km, 

x 2 = 7.9 km/s, Xi = x 3 = km/s (v cir ~ 6.31 km/s, t> par ~ 8.93 km/s). 

In an unperturbed case these values define an elliptic motion. 
Parameters of a potential are as follows: 

A_i = 0.04 km 4 /s 2 , Ai = 0.03 km 2 /s 2 , A 2 = -0.2 • 10~ 5 km/s 2 , 

B_i = 0.1 • 10" 4 km 4 /s 2 , B x = -0.0003 km 2 /s 2 , B 2 = 0.3 • 10~ 4 km/s 2 . 
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Figure 6: The case %a = 3, %b = 3. The hyperbolic type is in unperturbed case. 



Here the potential is not retaining. Coordinates of direction vector are b = (1, 1, 1) T . The 
roots of polynomials are as follows: 

6 « 2686, £a « 20699, Q\ « 4423 Q? G (6, 6), 

7^1^ 3256, i) 2)I)3 eC, Q°^ 5577 Q° > Vl . 
Therefore, the motion is unbounded. This is the case %a = 3, 13 = 4. The integration is 

0; 

-200000: 
Z -400000: 

-600000 ; 





Figure 7: The case %a = 3, = 4. 

carried during the time range corresponding approximately to t = 3.23 days. The particle 
trajectory is shown in fig. [71 

In the case of unbounded motion, to define the integration interval firstly one has to find 
the nearest pole of Q\{t) and/or Q%{t) in the direction of ascending r. Suppose this nearest 
pole is at r = T\. Then we choose a small positive value e and divide the segment [0,Ti — e\ 
into N equal subsegments. The value N is to be selected from practical reasons. The orbit 
should be visually smooth curve. In our examples the value iV = 100 was used. After that, 
the calculations by the formulae derived above are carried out in equidistant nodes. 

The following example demonstrates an application of our formulae for testing a numerical 
integration method. The original system of motion equations fl3"0l is considered. The Runge- 
Kutta-Fehlberg method of the eighth order with automatic choice of integration step is tested. 
The step is chosen by a method of the seventh order. The corresponding pair of programs, 
implemented in FORTRAN, is below noted as RKF8(7). Integration of equations (1301) 
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by RKF8(7) was performed with relative local error of the method e = 10~ 13 , and all 
calculations were carried out with double precision (real*8). The gravity parameter and the 
units of measurement are the same as above. A hypothetical particle is considered, repeatedly 
encountering the attracting centre. The trajectory obtained by explicit formulae is taken to 
be standard (reference). Its coordinates have been obtained using Maple with 32 digits (in 
FORTRAN this corresponds to quadruple precision (real*16)). 

Example 4. Initial values of coordinates and velocities of a particle are as follows: 

X\ = 7000 km, x 2 = km, £3 = 6000 km, 

x 2 = 7.9 km/s, X\ = x 3 = km/s {y cir ~ 6.58 km/s, v par ~ 9.30 km/s). 

In an unperturbed case we have an elliptic motion. 
Parameters of retaining potential are as follows: 

A_ x = 0.1 km 4 /s 2 , A 1 = -0.02 km 2 /s 2 , A 2 = -0.2 • 10" 5 km/s 2 

B_x = -0.004 km 4 /s 2 , B x = -0.001 km 2 /s 2 , B 2 = -0.001 km/s 2 

Coordinates of direction vector are b = (—1, —3, 1) T . The roots of polynomials $1 and $2 
are as follows: 

6~ 764, £3 ~ 58639, ^ 4459 Q? e (£2,6), 

V2 w 504, 773 w 7209, Q° 3 ^ 4761 Q° 3 G ( V 2, m)- 

Therefore, the motion is bounded. The case %a = 3, is = 3. 

The calculations were carried out during the time ranges corresponding to 1, 10, 50, 100, 
500, and 1000 revolutions of the particle around attracting centre without perturbations. 
The trajectory of the particle for three revolutions is shown in fig. [HI 




Figure 8: The case i& = 3, i% = 3. The motion is bounded. 

The table [1] contains the values, near the end of the integration interval, of the relative 
error for the energy constant 5H, the coordinates of particle position vector x«, and its 
absolute value r 
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Table 1: Estimation of the precision of numerical integration 



n 


t(day) 


5H x 10~ 12 


5x 1 x 10~ 12 


5x 2 x 10~ 12 


5x 3 x 10~ 12 


5r x 10" 12 


1 


.3382444 


1 


0.2 


1 


1 


0.4 


10 


4.9080991 


2 


6 


12 


213 


10 


50 


24.1940313 


41 


729 


104 


1667 


108 


100 


48.4322508 


53 


4399798 


523748 


95154 


330606 


500 


242.7821163 


294 


77898 


31418 


151259 


77206 


1000 


485.2955201 


556 


554500 


332688 


1067003 


330900 



where H° is the value at the initial moment, H at an arbitrary moment; x^ } r c are the values 
found by exact formulae. In the second column the intervals of physical time t (in days) are 
given, for which numerical integration of system fl30]) was carried out. 

From these data we can see that if the integration interval increases, the relative errors of 
H and 23 do not decrease. For coordinates X\, 22, and absolute value r, with n = 100, these 
errors increase, then they diminish, and then increase again. 

The numerical examples show efficiency of the formulae we obtained. Besides, the theorem 
3 allows to determine, given the initial position and velocity of a particle, whether its motion 
is bounded or unbounded. 

6. Conclusion 

In this paper we consider three sorts of coordinates (regular g-coordinates, bipolar coor- 
dinates, spherical coordinates). For each of the systems, the forms of potentials admitting 
complete separation of variables are given. Thus, the original equations for such potentials 
allow integration "in the sense of Sundman". In a similar way one can build, for regular q- 
coordinates, other coordinate systems for which Hamiltonian has orthogonal form, and with 
the use of Stackel theorem build potentials allowing the above-mentioned integrability. 

Application of these potentials is a separate and independent problem. Those potentials 
are of practical worth which approximate some real forces. 
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